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Summary 


Effects  of  large  and  small  scale  structure  on  the  propagation  of  Lg  and  regional  seismic 
phases  are  investigated  in  the  vicinity  of  Eurasian  seismic  arrays  KNET,  CNET,  and 
1LPA.  The  effects  of  small  scale  structure  are  predicted  using  models  of  Moho  and 
basement  topography  with  Lg  modeled  as  a  series  of  multiple  SmS  waves,  critically 
reflected  by  the  Moho.  Predictions  of  large  variations  in  Lg  amplitude  associated  with  the 
variations  in  crustal  thickness  surrounding  the  Hindu  Kush  and  Pamir  mountains  agree 
with  observations  at  KNET.  At  CNET  and  ILPA,  thick  basins  are  shown  to  have  large 
effects  on  shallow  crustal  sources,  blocking  either  a  portion  or  all  of  the  Lg  coda  for  paths 
crossing  these  basins.  The  effects  of  small  scale  structure  on  the  propagation  of -regional 
seismic  phases  are  investigated  using  a  Born  approximation  with  locked  mode  summation, 
with  scatterers  treated  as  equivalent  sources,  reradiating  the  incident  wavefield.  For 
shallow  sources,  it  is  shown  that  scattering  by  heterogeneity  close  to  the  source  and 
receiver  control  Lg  complexity,  equalizing  high  frequency  energy  on  all  components  of 
motion  regardless  of  source  type.  Experiments  with  data  from  the  Pinion  Flats  array 
suggest  that  anisotropic  distributions  of  heterogeneity  are  effective  in  simulating  the  coda  of 
regional  phases  in  regions  characterized  by  a  dominant  trend  of  deformation. 
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1  Introduction 


Whether  represented  by  critically  reflected  rays  at  the  Moho  (Gutenberg,  1955;  Bostock 
and  Kennett,  1990)  or  by  locked  modes  (Oliver  and  Ewing,  1957;  Knopoff,  Scwab  and 
Kausel,  1973;  Kennett,  1983;  Snieder,  1986;  Wang  and  Hermann,  1988;  Mitchell,  1995), 
the  energy  in  the  Lg  phase  is  the  usually  the  most  energetic  regional  seismic  phase  trapped 
in  the  crustal  waveguide,  and  hence  the  most  useful  for  nuclear  monitoring.  The  modal 
representation,  while  useful  in  1-D  models,  becomes  impractical  when  applied  to  2-D  and 
3-D  models  due  to  the  computational  expense  of  calculating  mode  coupling  coefficients. 
More  recently,  the  ray  approach  and  finite  difference  methods  are  used  in  the  study  of  Lg 
wave  propagation.  The  computational  expense  of  a  numerical  solution  to  the  wave  equation 
for  much  more  than  100  wavelengths  of  propagation,  however,  has  limited  finite 
difference  modeling  to  primarily  2-D  structure  and  low  frequency  bands. 

Ray  based  techniques  of  forward  modeling  are  sufficiently  fast  to  allow  predictions  of 
efficiency  of  the  propagation  of  Lg  and  are  well  suited  for  examining  the  attenuation  due  to 
the  scattering  and  blockage  of  Lg  in  the  perturbed  3-D  crustal  wave  guide.  Accurate 
application  of  the  ray  method  is  limited  by  the  ratio  A/s,  where  A  is  the  wavelength  and  s  is 
the  scale  length  of  the  medium.  Errors  in  the  prediction  of  body  wave  amplitudes  by  ray 
theory  are  typically  bounded  by  several  percent  when  the  non-dimensional  ratio  A/s  is  less 
than  several  0. 1  's.  The  smallest  scale  length  s  of  a  medium  is  either  the  minimum  of  the 
radius  of  curvature  of  a  discontinuity  or  the  minimum  of  the  quotient  of  velocity  (or 
density)  divided  by  the  magnitude  of  the  velocity  gradient  (or  density  gradient). 
Fortunately,  the  scale  lengths  of  structure  known  from  seismic  profiling  and  regional 
surveys  of  topography,  gravity  and  magnetic  fields  are  such  that  the  ratio  A/s  is  sufficiently 
small  in  at  least  a  portion  of  the  frequency  band  important  to  regional  monitoring  to  permit 
the  accurate  use  of  ray  theory  if  one  assumes  that  the  smaller,  unresolved,  scale  lengths  of 


structure  do  not  strongly  affect  the  propagation.  This  point  of  view  is  taken  in  the  first 
section  of  work  summarized  in  this  report  titled  Effects  of  Large  Scale 
Heterogeneity,  where  ray  theory  is  used  to  qualitatively  predict  paths  of  Lg  blockage  by 
large  scale  ( small  values  of  A/s )  variations  in  crustal  and  sedimentary  basin  thickness. 
Examples  of  qualitative  predictions  of  Lg  propagation  efficiency  are  compared  with  data 
observed  at  the  KNET,  CNET,  and  ILPA  arrays. 

As  the  ratio  A/s  becomes  much  larger  than  1 ,  paradoxically,  ray  theory  may  again  be  used 
accurately  if  scattering  by  small  scale  heterogeneities  is  approximated  by  the  Rayleigh-Born 
approximation.  If  the  velocity  and  density  perturbations  are  not  too  high,  a  heterogeneous 
medium  may  be  parameterized  by  a  sum  of  point  scatterers  and  the  effects  from  each 
scatterer  summed.  Since  the  ray  representation  of  Lg  is  already  given  as  a  sum  of  rays,  the 
requirement  of  an  additional  sum  over  ray  paths  to  and  from  scatterers  becomes  difficult  in 
both  computation  and  bookkeeping.  A  similar  type  of  problem  exists  in  the  conventional 
application  of  the  mode  representation  of  Lg  with  the  Bom  approximation,  where,  in 
theory,  every  mode  may  scatter  into  every  other  mode.  Although  many  of  these  mode 
interactions  can  be  shown  to  be  insignificant  in  energy,  the  approach  of  calculating  mode¬ 
coupling  coefficients,  while  useful  in  understanding  the  processes  of  scattering,  results  in 
an  inefficient  algorithm  for  waveform  synthesis  at  high  frequency.  In  the  second  section  of 
this  report  titled  Effects  of  Small  Scale  Heterogeneity,  an  alternative  algorithm  for 
the  mode-Bom  approximation  is  applied  to  small  scale  heterogeneity,  which  allows  it  to  be 
efficiently  applied  in  a  much  higher  frequency  band.  In  this  algorithm,  mode  coupling 
coefficients  are  not  computed.  The  complete  wavefield  is  calculated  at  each  point  scatterer 
and  reradiated  by  the  equivalent  source  representation  of  the  scatterer.  To  understand  the 
distribution  of  scatterers  responsible  for  equalizing  Lg  energy  on  three-components  of 
motion,  computational  experiments  are  made  with  various  depth  and  lateral  distributions  of 
scatterers  in  three-dimensions. 


2  Effects  of  Large  Scale  Structure 


2.1  Forward  modeling  method 

The  travel  times  and  amplitudes  of  multiple  SmS  rays  that  comprise  Lg  coda  are  computed 
by  integrating  kinematic  and  dynamic  ray  tracing  equations  (Cerveny  and  Hron,  1980; 
Cerveny,  1985).  The  Moho  and  the  boundary  between  a  sedimentary  layer  and  hard-rock 
layer  (basin  topography)  are  modeled  as  Epstein  transition  zones  (hyperbolic  tangents). 
With  the  Moho  and  basin  boundaries  so  modeled,  all  SmS  rays  become  turning  rays  in  the 
steep  gradient  region  of  the  Moho  transition  zone,  eliminating  the  need  for  computing  many 
reflection  coefficients  and  complex  descriptions  of  rays  interacting  with  three 
discontinuities  (surface,  basin/sediment  interface,  and  Moho).  The  effects  of  lateral 
variations  in  crustal  thickness  and  basin  thickness  on  Lg  efficiency  are  predicted  by  plotting 
SmS  ray  turning  and  bounce  points  or  cross  sections  of  the  crustal  waveguide  showing  ray 
paths.  Ray  bounce  points  may  be  connected  for  rays  shot  at  constant  vertical  and  varying 
azimuthal  take-off  angles  to  reveal  how  the  Lg  wavefront  propagates  at  distances  greater 
than  500  km.  Not  connecting  the  points  enables  closer  range  structure  effects  to  be 
observed.  The  existence  of  databases  of  Moho  and  basin  topography  sampled  on  grids  of 
10  km  or  smaller  spacing,  allows  this  approach  to  applied  to  predict  Lg  propagation 
efficiency  observed  at  several  arrays  in  Eurasia.  Results  are  reported  for  sources  within 
1000  km  of  the  KNET,  CNET  and  ILPA  arrays  (figure  1).  A  brief  geological  background 
of  these  regions  provides  a  context  in  which  to  study  the  regional  wave  propagation. 
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2.2  Tectonic  settings 


2.2.1  Garm  (KNET) 

The  Garm  region  (figure  2),  in  which  KNET  is  located,  is  one  of  the  most  seismically 
active  intra-continental  regions  in  the  world  and  is  characterized  by  intense  surface 
deformation.  Earthquake  focal  mechanism  studies  indicate  that  the  Garm  region  is 
dominated  by  nearly  horizontal  compression.  The  direction  of  principal  compression  varies 
from  approximately  north-west  in  the  western  region  to  nearly  north-south  in  the  North- 
East  region  (Lukk  and  Shevchenko,  1986,  1990;  Ramos,  1992).  Most  of  the  seismicity  is 
associated  with  the  active  deformation  of  the  sedimentary  rocks  in  the  Peter  the  First 
Range.  The  seismicity  may  be  associated  with  several  fault  zones,  the  most  prominent  is  a 
near  vertical  feature  beneath  the  axial  part  of  the  range  (Ramos,  1992).  Hypocenters  within 
the  Peter  the  First  Range  suggest  a  southeast  dipping  zone  of  seismicity  that  can  be  traced 
updip  to  the  Petrovsky  thrust  fault.  Under  this  upper  crustal  seismicity  at  a  depth  of  12  to 
17  km  is  an  aseismic  zone  that  may  be  related  to  low  strength  evaporites  at  the  base  of  the 
sedimentary  section.  The  aseismic  zone  is  then  underlain  by  an  anomalously  deep  crustal 
seismicity  at  a  depth  of  17  to  30  km.  This  deeper  seismicity  is  thought  to  be  associated  with 
a  southeast  dipping  zone  about  10  km  in  thickness,  of  intracontinental  subduction  whose 
surface  expression  is  the  Vakhsh  thrust  fault,  at  the  frontal  edge  of  the  Peter  the  First 
Range. 

Crustal  structure  in  the  Pamir-Tien  Shan  region  has  been  determined  by  inverting  for  a  1-D 
layered  velocity  model  using  data  from  two  networks  located  near  the  Soviet  earthquake 
prediction  laboratory  at  Garm,  Tadjikistan  (Hamburger  et  al.,  1993).  The  P  velocities  of 
two  upper  crustal  layers  are  indistinguishable  from  a  model  with  a  single  layer  over  a  half 


space.  Hambureger  et  ah’s  models  are  comprised  of  a  8  km  thick  upper  crustal  layer  of  5.2 
km/s  over  a  5.9  km/s  half  space,  and  a  5.2  km/s  layer  over  a  6.3  km/s  halfspace.  These 
results  contradict  previous  velocity  models  that  had  anomalously  low  velocity  zones  in  the 
mid-crust. 

2.2.2  Tibet  Plateau 

Figure  3  summarizes  the  broader  tectonic  setting  of  KNET.  Lg  has  been  observed  across 
most  of  Asia,  but  Lg  has  not  been  observed  for  paths  crossing  through  the  Tibetan  Plateau 
(Ruzaikin  et  al.,  1977).  This  has  been  attributed  to  several  factors,  the  change  in  crustal 
thickness  and/or  structural  discontinuities  at  the  boundaries  of  the  plateau  or  to  an  unusual 
velocity  structure  (absence  of  a  granitic  layer)  and  to  intrinsic  attenuation  within  the  plateau 
(Ruzaikin  et  ah,  1977;  Aki,  1980;  Gregerson,  1984;  Frankel  et  ah,  1990;  Mitchell,  1995). 
Crustal  thickness  in  the  Tibetan  Plateau  is  greater  than  the  surrounding  tectonic  regions  of 
the  Indian  Shield  to  the  south  and  the  Tarim  Basin  to  the  north  (Molnar,  1988).  This 
significantly  deforms  the  crustal  waveguide.  Complicated  structures  of  the  Himalayan 
boundary  thrust  and  the  northern  Kunlun  fault  may  also  cause  significant  scattering  of  Lg 
energy. 

The  Tibetan  Plateau  is  bounded  by  the  Kunlun  Mountains  in  the  north  and  the  Himalayas  to 
the  south.  It  has  an  average  elevation  of  5  km  and  an  areal  extent  of  nearly  700,000  square 
kilometers.  The  uplift  of  the  plateau  is  the  result  of  the  collision  of  the  Indian  plate  with 
Eurasia  that  began  in  the  Middle  Eocene  approximately  45  my  (Molnar  and  Tapponier, 
1975).  Relative  motion  of  the  plates  continues  at  a  rate  of  5  cm/yr.  The  Tibetan  Plateau  is 
highly  fractured  due  to  the  suture  zones  bounding  the  many  terranes  that  make  up  the 
plateau  and  the  abundant  normal  and  strike-slip  faults  in  the  uppermost  20  km  along  which 
the  crust  is  actively  deforming  (Randall  et  al.,  1995).  Most  earthquakes  that  occur  within 


the  plateau  are  located  in"  the  top  5-10  km  of  the  crust  (Molnar  and  Chen,  1983;  Zhao  and 
Helmberger,  1991).  Geophysical  observations  of  shallow  seismicity  are  consistent  with 
anomalously  high  heat  production  beneath  the  northern  portion  of  the  plateau.  The 
earthquake  cutoff  temperature  is  estimated  to  be  between  250  and  450  C  (Chen  and  Molnar, 
1983).  Volcanic  flows  of  both  granitic  and  basaltic  composition  are  observed  at  the  surface 
throughout  the  northern  section  of  the  plateau  (Dewey  et  al„  1988).  This  indicates  a  mantle 
source  of  volcanism  causing  crustal  heating  as  it  propagates  to  the  surface.  Crustal  heating 
is  likely  to  significantly  increase  attenuation  (decrease  Q)  (Frankel  et  al.,  1990).  High 
crustal  temperatures  are  not  clearly  understood  and  may  be  due  to  the  concentration  of 
radiogenic  crustal  material  resulting  from  the  crustal  shortening  during  the  plateau 
formation  (Zhao  and  Helmberger,  1991)  or  due  to  a  high  heat  flux  from  the  upper  mantle 
(Molnar,  1988). 

2.2.3  Caucasus  and  Iran  (CNET  and  ILPA) 

The  Caucasus  region,  in  which  CNET  (figures  4a-c  )  is  located,  is  a  mixture  of  tectonic 
units  created  by  the  collision  of  the  Arabian  and  Eurasian  plates  that  initiated  approximately 
25  million  years  ago.  The  interaction  of  these  large  plates  determines  the  meridional 
direction  of  maximum  compression  stresses  in  the  Caucasus  and  kinematics  of  movements 
in  the  system  of  large  and  small  faults,  e.g.,  North  Anatolian,  East  Anatolian,  Zagros, 
Pambak-Sevan,  etc.  The  Caucasus  orogeny  is  similar  to,  but  younger  than  the  Himalayan 
orogeny,  with  the  orientation  of  the  colliding  Arabian  plate  the  same  as  the  north-east 
verging  Indian  plate).  The  Lesser  Caucasus  (41.5  latitude)  consists  of  a  north-verging 
system  of  folds  that  have  their  origin  in  the  continental  collision.  The  Greater  Caucasus 
(42.5  latitude)  was  formed  from  south-verging  system  of  folds,  related  to  the  northward 
thrust  of  the  Arabian  plate.  These  different  tectonic  units  are  separated  by  the 
Transcaucasian  Massif  (Adamia  et  al.,  1991).  The  active  tectonic  setting  is  a  compressive 


one  with  evident  slip  fronts  and  folded  structures.  Folds  are  locally  interrupted  by 
continuous  structural  lineations  that  suggest  a  strike-slip  movement.  The  most  significant 
lineament  is  parallel  to  the  orogenic  belt. 


The  Greater  Caucasus  has  three  large  morphostructural  units:  the  Kura  basin,  the  Greater 
Caucasus,  and  the  foothills.  The  deformations  of  the  Greater  Caucasus  involve  the  whole 
crust  and,  most  likely,  the  upper  mantle,  all  the  way  to  the  Lesser  Caucasus.  The  general 
trend  of  large  morphostructures  from  west  to  north-west  is  oblique  to  the  direction  of  the 
prevailing  compression.  The  Greater  Caucasus  region  is  highly  faulted  and  deformed  due 
to  mountain  building.  The  meridional  faults  faults  of  the  Trans-Caucasian  extension  are 
extensional  faults,  the  near-latitudinal  Caucasian  faults  are  compression  upthrow  shear 
thrust  faults  and  upthrow  shifts,  and  the  diagonal  faults  extended  north-west  are  shear 
structures,  i.e.,  strike-slip  displacements  (Rogozhin  et  al.,  1993). 

The  Lg  phase  propagates  well  throughout  most  of  the  region  except  for  events  coming  from 
the  northern  Black  Sea  and  the  Caspian  Sea  regions  and  events  in  the  Hormuz  region  of 
Iran  (at  the  eastern  end  of  the  Persian  Gulf  opposite  Oman)  (Sweeney,  1996).  The  Black 
Sea-Caspian  Sea  Lg  blockage  is  speculated  to  be  related  to  a  change  in  crustal  thickness 
beneath  the  Black  Sea  and  the  Caspian  Sea  or  to  basin  effects,  or  a  combination  of  both. 
Sediments  are  as  thick  as  20  km  beneath  the  Caspian  sea.  Blockage  of  Lg  from  the 
Hormuz  region  is  not  understood,  as  Lg  is  not  blocked  from  other  similar  parts  of  the 
Zagros  range  to  the  northwest.  The  Kazerun  line  is  a  N-S  trending  line  consisting  of  the 
Borazjan,  Kazerun,  and  High  Zagros  faults  and  divides  the  Zagros  range  (a  fold-thrust  belt 
lying  to  the  southwest  of  the  Main  Zagros  Reverse  Fault)  into  two  different  regions.  In  the 
southeast  there  is  a  thick  ( 1  km)  late  Precambrian  evaporite  unit  that  has  been  mobilized  into 
extensive  salt  diapirs  and  nappes;  to  the  northwest,  the  salt  unit  is  much  less  prominent  and 


mobilized  salt  is  much  rarer.  The  mobile  salt  is  thought  to  be  a  possible  cause  of  scattering 
of  Lg  in  this  region.  In  the  extreme  southeastern  comer  of  Iran  there  is  oceanic  subduction 
and  local  magma  formation,  with  an  inboard  volcanic  arc  delineated  by  a  string  of 
Quaternary  volcanic  centers.  Quaternary  volcanism  is  also  associated  with  the  central 
Iranian  plateau  southwest  of  the  city  of  Meshad. 

2.3  Instrumentation  and  data 

2.3.1  KNET 

The  Kyrgyz  Seismic  Network  consists  of  10  stations  and  is  located  near  the  capital  of 
Kygyzstan,  Bishkek  (figure  5).  Tectonically,  it  is  located  approximately  150  km  north  of 
the  Garm  region  and  the  western  flank  of  the  Himalayan  plateau.  Each  station  is  comprised 
of  three  component  Streckeisen  STS-2  seismometers  with  Reftek  RT72A-02  dataloggers 
with  100  sps  and  20  sps  datastreams  for  three  high  gain  and  three  low  gain  channels.  The 
KNET  stations  are  all  located  in  hard  rock,  but  stations  CHM,  TKM,  KBK,  and  EKS2 
have  shallow  sediments  located  within  100  meters  of  their  vaults. 


2.3.2  CNET 

In  August  1990,  a  Lamont-IRIS  team  deployed  a  seismic  network  in  the  Caucasus.  The 
original  intent  was  to  span  the  seismically  active  Caucasus,  but  the  array  size  was  reduced 
due  to  political  turmoil  in  the  region.  The  scope  was  limited  to  a  small  aperture  array  in  the 
region  near  Kislovodsk,  on  the  northern  slopes  of  the  Greater  Caucasus  and  on  the 
southern  Russian  Platform  (figure  4). 
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2.3.3  ILPA 


The  set  of  continuous  waveform  data  from  the  Iranian  Long-Period  Array,  (ILPA  ,  figure 
6)  used  in  this  study  was  compiled  by  Jerry  Sweeney  (1996).  The  ILPA  array  was 
operational  during  the  late  70’s.  The  “Long  Period”  array  contained  both  short  period  and 
long  period  instruments.  The  array  sensor  sites,  remotely  located  with  respect  to  the 
recording  station,  are  approximately  60  kilometers  southwest  of  Tehran,  Iran.  Stations  2 
through  7  are  arranged  in  a  hexagon  configuration  with  Station  1  located  in  the  center  of  the 
array.  The  sites  are  all  located  in  rocky  areas  where  very  little  overburden  exists,  and  on 
high  points  with  respect  to  the  surrounding  area. 

2.4  Results  of  ray  experiments 

2.4.1  KNET 

Predictions  of  Lg  propagation  efficiency  computed  by  ray  theory  are  compared  with  the 
waveforms  of  specific  events  chosen  to  illustrate  the  effects  of  varying  Moho  topography. 
Of  particular  interest  are  the  effects  of  Moho  topography  on  Lg  paths  that  cross  extreme 
lateral  gradients  in  crustal  thickness  associated  with  the  Pamir  Range  and  the  Himalayan 
Plateau  as  well  as  Lg  paths  that  cross  the  relatively  flat  Tien  Shan  region  to  the  east  of 
KNET. 

A  one  dimensional  velocity  model  (JSPC,  1993)  for  the  region  is  used  along  with  the 
database  of  Moho  topography  by  Fielding  et  al.  (1992)  to  provide  a  three  dimensional 
model  for  the  dynamic  ray  tracing.  The  variations  in  Lg  coda  from  regional  seismograms 
recorded  by  the  KNET  array  in  the  Republic  of  Kyrgyzstan  are  compared  with  predicted 
variations  in  the  paths  of  multiple  SmS  waves  comprising  the  Lg  phase.  Ray  trajectories, 
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SmS  turning  points  and  surface  bounce  points  (Bostock  and  Kennett,  1990),  plotted  in 
map  view,  are  used  in  the  interpretation  and  prediction  of  structural  effects  on  Lg. 

Events  were  selected  based  on  the  following  criteria:  shallow  depth,  a  good  signal  to  noise 
ratio,  range  from  stations  between  300  and  1200  km,  and  source-receiver  paths  that  cut 
through  regions  where  the  Moho  depth  varies  greatly.  The  data  are  processed  using  the 
Seismic  Analysis  Code  (Lawrence  Livermore  National  Laboratory,  1990).  The  data  are 
band  passed  (0.5-5  Hz)  using  a  Butterworth  filter  with  four  poles  and  then  plotted  in  a 
record  section  for  each  event.  Approximately  40  events  met  the  selection  criteria.  Out  of 
that  data,  three  events  were  chosen  (figure  7a,b)  that  represented  sources  from  each  of  the 
three  tectonic  regions  near  KNET.  Event  1  is  located  in  the  Pamir  region.  Event  2  is  in  the 
Himalayan  Plateau,  and  event  3  is  located  in  the  Tien  Shan  region  to  the  east. 


Event  1.  6  May  92.  This  event  was  chosen  for  the  ray  paths  that  pass  through  a  region 
north  of  the  Pamir  Mountains  where  the  depth  of  the  Moho  increases  to  70  km  and  is 
characterized  by  steep  gradients.  The  array  stations  had  a  back  azimuth  range  of  198.02  to 
209.86  degrees  and  ranged  in  distance  from  527.9  to  61 1 .3  km.  The  effects  of  the  steep 
Moho  gradient  are  seen  in  the  Lg  amplitude  and  duration,  which  differ  greatly  across  the 
array  (figure  8).  Stations  TKM  (baz  209.86;  dist.  527  km)  and  CHM  (baz  204.86;  dist 
592  km)  have  very  little  Lg,  while  station  AAK  (baz  204.47;  dist  546  km)  has  noticeable 
Lg  and  EKS2  (baz  198.02;  dist  527  km)  has  significant  Lg.  Ray  tracing  predicts  strong 
variations  for  rays  passing  this  region.  The  map  view  of  SmS  ray  bounce  points  (figure  9) 
for  this  source  receiver  pair  shows  that  the  wave  front  is  distorted  by  the  second  bounce 
and  a  two  dimensional  caustic  develops  that  transects  the  array.  This  correlates  well  with 
the  differences  seen  in  Lg  from  one  side  of  the  array  to  the  other.  Two-dimensional  vertical 
cross  sections  of  ray  paths  projected  onto  vertical  planes  through  back  azimuths  for  this 
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event  (figure  10)  show  significant  effects  of  the  rapid  lateral  variation  in  crustal  thickness 
The  ray  trajectories,  which  differ  only  slightly  in  azimuth,  differ  dramatically  in  focusing 
and  defocusing  at  distances  as  close  as  150  km.  The  rays  diverge  after  the  first  bounce  and 
the  range  of  coherency  of  the  Lg  wavefront  is  much  less  than  150  km. 

Event  2.  27  June  92.  Stations  in  the  array  had  a  back  azimuthal  range  of  139.17  to  146.57 
degrees  and  ranged  in  distance  from  987  to  1047  km.  The  epicenter  for  this  event  is  in  the 
thickened  crust  (70  km)  associated  with  the  Hindu  Kush  region.  The  record  section  (figure 
11)  shows  strong  variations  across  the  array  in  Lg  coda  shape,  duration  and  amplitude  ratio 
to  Pg,  which  is  expected  for  paths  that  cross  regions  of  strong  Moho  topography.  The  Lg 
coda  at  station  EKS2  has  a  different  shape  and  a  lower  Lg/Pg  ratio  than  the  rest  of  the 
array.  The  ray  tracing  for  this  event  shows  most  of  the  array  lying  in  a  region  of 
defocusing  (figure  12),  while  the  western  most  station  (EKS2)  lies  on  a  two  dimensional 
caustic  or  region  of  apparent  focusing.  This  may  account  for  the  fuller  Lg  coda  seen  at 
station  EKS2.  The  cross  section  of  ray  paths  through  the  crust  for  this  event  (figure  13) 
illustrates  the  effects  of  focusing  and  defocusing  of  rays  along  this  path. 

Event  3.  25  Mar  92.  This  event  traverses  a  region  of  Moho  east  of  the  array  that  is 
relatively  flat  and  is  approximately  47  km  deep.  The  record  sections  (figure  14a)  show 
similarities  in  both  Lg  coda  length,  shape  and  amplitude  ratio  to  Pg.  The  similarities  in  Lg 
coda  are  consistent  with  the  expected  effect  of  the  flat  Moho  between  the  source  and  the 
array.  Ray  tracing  shows  that  the  rays  passing  through  this  region  are  not  greatly  affected 
by  the  slight  changes  in  Moho  topography.  The  ray  tracing  results  are  shown  in  map  view 
of  SmS  bounce  points  (figure  1 5)  and  two  dimensional  cross  sections  of  ray  paths  (figure 
16).  The  bounce  points  show  a  relatively  coherent  wavefront  reaching  the  array.  This  is 
also  illustrated  in  the  cross  section  of  ray  paths.  The  cross  section  of  ray  paths  for  this 


event  show  the  even  spacing  of  rays  that  one  expects  for  rays  bouncing  between  nearly 
parallel  surfaces.  The  Lg  rays  in  this  case  remain  coherent  out  to  a  range  of  600  km. 

The  ray  synthetics  (figures  6b,  1  lb,  14b)  are  calculated  for  events  1-3.  The  ray  synthetics 
show  large  variation  for  azimuths  that  differ  by  as  little  as  one  degree  in  regions  of  steep 
Moho  gradients.  The  synthetics  for  events  1  and  2  differ  in  content  and  amplitude  as  a 
result  of  ray  tracing  through  Moho  topography  associated  with  the  Hindu  Kush  and  Pamir 
mountain  ranges.  The  synthetics  simulating  event  3  have  similar  amplitudes  and  arrivals 
across  the  array  as  a  result  of  ray  tracing  through  relatively  flat  Moho  topography 
associated  with  the  region  east  of  the  array.  Although  there  is  little  hope  with  a  simple 
homogeneous  crust  of  reproducing  the  coda  complexity  seen  in  the  data,  the  time  spacing 
and  amplitudes  of  individual  SmS  rays  are  useful  in  assessing  affects  on  the  complete  Lg 
coda. 

Reflectivity  (Fuchs  and  Muller,  1971)  synthetics  are  used  to  model  event  #3  from  the  Tien 
Shan  region  to  the  east,  where  there  is  little  variation  in  crustal  thickness.  A  velocity  model 
is  assumed  from  the  JSPC  technical  note,  with  Moho  depth  taken  from  mean  values  along 
this  path  in  the  Fielding  database.  Using  this  structural  model,  the  Lg  to  Pg  amplitudes  do 
not  match  well.  A  simple  trial  and  error  inversion  is  performed  by  comparing  the  synthetics 
varying  Q  values  to  fit  the  waveforms  from  event#3  (figure  17).  .  The  results  from  a 
preliminary  crustal  model  and  Q  values  similar  to  those  in  New  England  (Pulli,  1984) 
show  little  resemblance  to  the  event  #3  and  suggest  that  the  event  may  be  mislocated. 
Synthetics  assuming  a  location  100  km  closer  to  the  array  provide  a  better  match  to  the 
observed  waveform,  suggesting  that  mislocation  can  be  a  severe  impediment  to 
interpretation  of  propagation  from  some  of  the  events  in  this  region. 
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2.4.2  CNET 


While  specific  paths  may  have  particular  interest,  a  broader  assessment  of  regional 
propagation  can  be  obtained  by  shooting  the  rays  outward  from  a  station.  In  this  way  it  is 
possible  to  make  propagation  predictions  over  the  entire  360  degree  azimuthal  range  from 
an  array,  including  paths  from  regions  in  which  natural  seismicity  is  absent.  The 
justification  for  doing  this  lies  in  the  reciprocity  theorem  for  the  elastic  Green’s  function. 
Figure  18  shows  the  triggered  events  within  10  degrees  of  CNET  for  1992.  Regional  data 
observed  at  CNET  has  an  azimuthal  dependency.  Events  located  to  the  east  of  CNET  are 
observed  to  exhibit  strong  Lg  excitation.  Figure  19a  shows  a  typical  event  with  strong  Lg 
from  the  region  to  the  east.  Events  to  the  southwest  of  CNET  are  observed  to  have 
relatively  weak  or  unidentifiable  Lg,  this  has  also  been  noted  in  a  study  of  Lg/P  ratios  by 
Abers  (1995).  Figure  19b  shows  a  typical  event  from  the  south-southwest  that  exhibits 
weak  or  absent  Lg.  Crustal  and  basin  thickness  variations  in  the  vicinity  of  CNET  and 
ILPA  taken  from  the  Fielding  data  set  are  illustrated  in  figure  20.  CNET  lies  near  a  strong 
gradient  in  Moho  thickness  associated  with  the  Caucasus  and  other  strong  gradients 
associated  with  the  Black  Sea  to  the  southwest. 

The  ray  study  in  this  region  is  centered  on  station  KIV  in  CNET.  Predictions  of  Lg 
efficiency  are  sought  by  modeling  the  region  using  variations  in  crustal  thickness  as  well  as 
variations  in  basin  thickness.  Rays  are  shot  at  dip  angles  of  60,  70  and  80  degrees. 
Sources  are  located  at  both  shallow  (1  km)  and  mid  crustal  (15  km)  depths.  A  map  view  of 
bounce  points  for  the  model  that  has  only  crustal  thickness  variations  is  plotted  in  figures 
21a,c-23a,c.  The  Fielding  data  set  for  basin  thickness  was  then  implemented  in  the  velocity 
model  and  the  same  tests  were  run.  A  map  view  of  bounce  points  for  rays  in  the  CNET 
model  with  both  crustal  and  basin  thickness  variations  is  plotted  in  figures  21b,d-23b,d. 


The  large  scale  effects  of  the  three-dimensional  varying  surfaces  of  the  Moho  and  basins  on 
Lg  efficiency  are  strongly  correlated  with  steep  gradients  in  those  surfaces.  In  regions  of 
steep  Moho  gradients  only  (no  basins),  rays  are  bent  and  blocked  significantly,  and  the 
trends  in  Lg  efficiency  are  fairly  easy  to  follow  with  the  eye.  Strong  defocusing  is  predicted 
to  the  west  and  southwest  of  CNET,  where  events  showing  weak  or  unidentifiable  Lg  are 
located.  Adding  basins  to  the  ray  modeling  greatly  amplifies  the  effects  from  the  Moho  and 
decreases  the  coherency  of  the  Lg  wavefronts.  Some  of  rays  responsible  for  segments  of 
Lg  coda  may  be  trapped  in  the  basins  (figures  21b,d-23b,d),  where  there  may  be  local 
amplification.  For  shallow  events  of  1  km  depth  (21b-23b),  basins  play  a  major  role  in 
scattering  and  blocking  the  rays,  as  the  rays  start  out  trapped  in  the  basins  and  may  get 
blocked  at  the  basin  crust  boundary.  For  the  deeper  events  of  15  km,  Lg  rays  in  the  model 
having  only  Moho  variations  seem  to  propagate  more  coherently  than  the  those  from  the 
shallow  crust.  Deeper  events  in  the  model  with  both  basin  and  Moho  produce  rays  that  may 
propagate  100  km  before  getting  trapped  in  the  basins.  For  many  azimuths,  the  Lg  does  not 
get  trapped  in  the  basins  at  all. 

For  the  CNET  array,  the  effects  due  to  the  thinning  of  the  crust  under  the  Black  Sea  is 
accentuated  by  thick  sediments  and  rays  are  strongly  scattered  and  blocked.  This  is  seen  in 
the  map  view  plots  of  both  ray  paths  and  of  SmS  bounce  points  (figures  21-23).  The  ray 
path  plots  illustrate  the  three  dimensional  nature  of  the  scattering,  where  some  rays  undergo 
severe  bending  over  small  distances  of  100  to  200  km. 

2.4.3  ILPA 

The  observations  for  Lg  used  in  this  study  are  those  of  Sweeney  (1996).  Sweeney, 
reviewed  the  200  plus  events  in  the  data  set,  performing  seismic  phase  picking  and  data 


quality  analysis.  Phases  were  identified  as  present  only  if  the  phase  was  distinct  and 
unambiguous.  Phases  identified  as  Lg  were  timed  on  the  vertical  channels.  Sweeney’s  plot 
of  the  Lg/P  ratios  are  shown  in  figure  6.  In  the  figure,  events  for  which  a  phase  is  present 
are  plotted  using  a  filled  dot;  where  the  indicated  phase  is  not  present,  the  event  location  is 
plotted  using  a  doughnut.  The  propagation  path  in  each  case  is  from  the  event  to  the 
location  of  the  array.  The  locations  of  important  faults  and  the  locations  of  Quaternary 
volcanic  centers  (triangles)  are  noted  as  well. 

Rays  used  in  the  ILPA  were  calculated  as  in  the  CNET  array.  A  source  was  located  at  a 
single  station  and  rays  were  shot  for  depths  of  1  and  15  km,  as  well  as  for  take  off  angles 
of  60  ,70,  and  80  degrees.  A  map  view  of  bounce  points  for  rays  in  the  ILPA  model  with 
both  crustal  and  basin  thickness  variations  as  well  as  variations  in  source  depth  and  take  off 
angle  is  plotted  in  figures  24-26. 

The  Moho  (without  basins)  effect  is  moderate  with  blockage  and  scattering  to  the  north  and 
north-west:  (figure24a,c-26a,c).  Propagation  past  500  km  is  generally  highly  incoherent. 
Near  ILPA  the  Moho  is  relatively  flat  and  it  is  the  basins  that  do  the  most  scattering  of  Lg. 

A  main  feature  that  blocks  rays  to  the  north  of  the  array  for  both  the  shallow  and  deep 
sources  is  the  deep  basin  that  is  beneath  the  Caspian  Sea.  This  ray  blockage  is  seen  from 
38E  to  42N  to  55E  to  38N  (figures  24-26)  and  agrees  well  with  the  low  Lg  to  P  ratios  seen 
by  Sweeney  (1996)  (figure  24)  and  with  prior  studies  of  Wallace  et  al,  (1992)  who 
observed  inefficient  Lg  propagation  across  northern  Iran  and  the  southern  Caspian  region, 
and  with  Kadinsky-Cade  et  al.  (1981)  who  found  that  for  travel  paths  across  the  southern 
Caspian  region,  the  Lg  phase  is  scattered  or  blocked.  Maps  of  bounce  points  from  models 
with  basins  differ  greatly  between  a  shallow  and  mid-crustal  source.  This  is  due  to  the 
shallow  sediments  to  the  south  and  east  of  the  array  and  may  provide  a  means  of  depth 
discrimination  in  the  area.  The  northwest  trend  in  trapped  bounce  points  seen  south  of  the 
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array  from  the  models  with  basins  included  (figures  24b,d-26b,d)  may  be  associated  with 
the  plate  boundary  and  the  Zagros  fault.  This  correlates  well  with  the  low  Lg/Pg  ratios  seen 
in  Sweeney’s  plot  south  of  the  ILPA  array. 

2.5  Summary  of  ray  experiments  on  Lg  propagation 

2.5.1  Crustal  thickness  variations 

The  most  important  effects  of  large  scale  structure,  including  strong  Lg  attenuation  and 
blockage,  can  be  estimated  from  deterministic  crustal  models  based  on  regional  refraction 
and  deep  seismic  reflection  surveys,  surficial  topography,  and  gravity.  Crustal  and  basin 
thickness  variations  have  been  incorporated  into  ray  tracing  models  for  the  purposes  of 
modeling  regional  Lg  wave  propagation.  Ray  tracing  predicts  the  strong  regional  variations 
seen  in  Lg  efficiency  due  to  strong  Moho  gradients  seen  in  data  collected  at  the  Kyrgyz 
array.  The  plots  of  the  SmS  bounce  points  for  the  events  in  this  study  suggest  that  a  large 
transition  in  Lg  amplitudes  occurs  for  paths  traversing  regions  of  strong  Moho  gradient 
associated  with  the  Hindu  Kush  and  Pamir  mountains.  An  effect  of  crustal  thickness 
variation  is  a  modification  of  SmS  ray  paths  that  results  in  either  a  compression  or 
stretching  of  Lg  coda.  In  regions  of  uniform  crustal  thickness,  however,  reflectivity  CRFL 
synthetics  yield  reasonably  good  fits  to  observed  data.  For  example,  in  the  Tien  Shan 
region  east  of  the  KNET  array,  the  Moho  topography  is  reasonably  flat,  and  one¬ 
dimensional  modeling  yields  new  estimates  for  seismic  Q  that  are  slightly  lower  than  the 
previous  values  for  the  region. 


2.5.2  Effects  of  basins 


Incorporation  of  sedimentary  basins  in  ray  modeling  of  Lg  at  CNET  and  ELPA  shows  that 
thick  basins  may  have  large  effects  of  shallow  crustal  sources,  blocking  a  either  a  portion 
or  all  of  the  Lg  phase  along  paths  crossing  these  basins.  The  basin  effects  are  strongly 
sensitive  to  source  depth,  nearly  vanishing  for  sources  in  the  deep  crust.  In  some  regions, 
the  depth  sensitivity  of  the  basin  effect  may  provide  the  basis  for  the  construction  of  an 
event  discrimination  by  source  depth  from  Lg  amplitudes. 
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Figure  1 :  This  topographical  map  of  Eurasia  illustrates  the  location  of  the  three  seismic 
arrays  in  the  regions  studied  (ETOP05,  1988). 
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Figure  2:  a)  Regional  geology  and  b)  seismicity  and  major  faults  of  the  Garm  region.  (M. 
Hamburger,  Geophys.  J.  Int.,  1993) 
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Figure  3:  Map  of  the  Himalaya  and  central  Asia,  based  on  Searle’s  (1996)  interpretation 
of  Avouac  and  Tapponnier  (1993)  showing  zones  of  presently  active  crustal  shortening 
and  estimated  kinematic  rates  of  shortening,  extension,  and  strike-slip  fault  motion. 


Latitude 


21 


Figure  4a:  Topography  in  the  vicinity  of  the  CNET  and  ILPA  arrays,  with  an  exploded 
view  (JSPC  1992)  of  the  Caucasus  array. 
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Figure  4b:  CNET/Caucasus  faults  and  tectonics. (Rhogozhin  and  Phillip,  1991) 
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Figure  4c:  Tectonic  setting  and  major  active  faults  near  CNET  and  ILPA  arrays  (Spitak 
Earthquake  database  for  1997) 
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Station  Code 

Station  Name 

Latitude  (deg) 

Longitude  (deg) 

—Mil 

AAK 

Ala-Archa 

42.6333 

74.4944 

1.6800 

AML 

Almalyashu 

42.1311 

73.6941 

3.4000 

BGK2 

Belogorka 

42.6451 

74.2274 

1.6400 

CHM 

Chumysh 

42.9986 

74.7513 

0.6550 

EKS2 

Erkin-Sai 

42.6615 

73.7772 

1.3600 

KBK 

Karagaibulak 

42.6564 

74.9478 

1.7600 

KZA 

Kyzart 

42.0778 

75.2496 

3.5200 

TKM 

Tokmak 

42.8601 

75.3184 

0.9600 

UCH 

Uchtor 

42.2275 

74.5134 

3.8500 

USP 

Uspenovka 

43.2669 

74.4997 

0.7400 

Table  1.  Kyrgyz  stations  and  locations  (JSPC,  1993). 


Figure  5:  Kyrgyz  seismic  telemety  network  (JSPC,  1992). 


24 


A  location  of  Quaternary  volcanic  center 
•  location  of  event  with  L  g  phase  present  at  ILPA 

O  location  of  event  with  L  g  phase  not  present  at  ILPA 

yzRF  Main  Zagros  reverse  fault  B  Borazjan  fault  k  Kazerun  fault 
hzf  High  Zagros  fault  (faults  after  Pig.  1  of  Berberian,  1  995) 
mdf  Makran  Reformation  front  (oceanic  subduciion*  after  Byrne  at  ah,  1992) 


Figure  6:  Interpretation  of  ILPA  data  (1978-1979):  Lg  phase  propogation  and  location  of 
ILPA  and  Lg  blockage.  Note  the  Zagros  fault  system.  (J.  Sweeney ,  1996) 
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Moho  topography  (crustal  thickness)  based  on  (crustal  thickness)(Fielding,  1992)  with  3 
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Figure  9:  Map  view  of  ray  bounce  points  from  dynamic  ray  tracing  for  event  #1.  The 
wavefront  is  strongly  affected  by  the  steep  Moho  topography.  A  two  dimensional  caustic 
appears  to  be  forming  down  the  middle  of  the  array. 
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Figure  11a:  Record  section  of  plot  event  #2  in  figure  Figure  1  lb:  Ray  synthetic  of  event  2  to  stations  AAK,  CHM,  and 

7a,  27  June  92.  Data  was  bandpassed  filtered  1-5  Hz.  TKM.  The  synthetics  for  the  different  stations  show  the  effects  of 

The  Lg  coda  varies  across  the  array.  the  Moho  topography  associated  with  the  Hindu  Kush. 
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Figure  13:  Rays  for  event  #2  at  take  off  angles  63-71  degrees  for  paths  to  stations  TKM 
and  EKS2.  Both  paths  show  development  of  unique  focusing  and  defocusing  of  rays  due 
to  the  gradients  in  the  Moho  topography. 
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Figure  16:  Tien  Shan  Vertical  Cross-section.  Rays  for  event  #3  at  take  off  angles  61-69 
degrees  for  paths  to  stations  AAK  and  CHM.  Due  to  the  relatively  flat  Moho  along  the 
ray  path,  ray  paths  are  similar  out  to  500-600  km. 


Figure  17a:  Seismograms  from  event  #3  recorded  at  stations  AAK  and  CHM.  Band 
passed  l-5Hz. 


Figure  17c:  At  ranges  of  420-430  km  there  is  a  good  match  in  coda  length,  shape  and 
Lg/Pn  amplitude,  indicating  a  need  for  relocation. 
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All  events  plotted  40°  45°  50° 

EH  Event  in  the  PDE  T 

■  Event  not  in  the  PDE  Longitude 


#  Ho  Mb 

■  Mb  =  2.0 

*  Mb  =  3.0 

■  Mb  =  4.0 

■  Mb  =  5.0 


Figure  18:  Events  triggered  during  January  1,  1992  to  November  9,  1992  from  CNET, 
centered  at  KTV,  circles  indicate  increments  of  5°  (JSPC,  1994). 
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Figure  19a:  Example  of  efficient  Lg  propagation  from  an  event  to  the  east  of  KNET. 
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Figure  19b:  Example  of  inefficient  Lg  propagation  from  an  event  to  the  west  and  south 
west  of  KNET. 


Depth  in  kilometers 

Figure  20:  Moho  structure  in  the  vicinity  of  CNET  and  ILPA  arrays  based  on  Fielding  et 
al„  1993. 
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Figure  21a:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  database.  Ray  take-off  dip  angle  is  60°,  source  depth  is  1  km. 

Rays  were  shot  at  1°  azimuthal  increments  for  360°. 


Figure  21b:  Map  view  of  ray  boi 
Fielding  Moho  and  basin  databa 
Rays  were  shot  at  1°  azimuthal  i 
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Figure  21c:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  database.  Ray  take-off  dip  angle  is  60°,  source  depth  is  15km. 

Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  2 Id:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  and  basin  databases.  Ray  take-off  dip  angle  is  60°,  source  depth  is  15  km 
Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  22a:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  database.  Ray  take-off  dip  angle  is  70°,  source  depth  is  1  km. 

Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  22b:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  and  basin  databases.  Ray  take-off  dip  angle  is  70°,  source  depth  is  1  km. 
Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  22c:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  database.  Ray  take-off  dip  angle  is  70°,  source  depth  is  15  km. 
Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  22d:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  and  basin  databases.  Ray  take-off  dip  angle  is  70°,  source  depth  is  15  km. 
Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  23a:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  database.  Ray  take-off  dip  angle  is  80°,  source  depth  is  1  km. 

Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  23b:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  and  basin  databases.  Ray  take-off  dip  angle  is  80°,  source  depth  is  1  km. 
Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  23c:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  database.  Ray  take-off  dip  angle  is  80°,  source  depth  is  15  km. 

Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  23d:  Map  view  of  ray  bounce  points  from  a  regional  model  incorporated  with 
Fielding  Moho  and  basin  databases.  Ray  take-off  dip  angle  80°,  source  depth  15km. 
Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  24a:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  database.  Ray  take-off  dip  angle  is  60°,  source  depth  is  1  km. 

Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  24b:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  and  basin  databases.  Ray  take-off  dip  angle  is  60°,  source  depth  is  1  km 
Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  24c:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  database.  Ray  take-off  dip  angle  is  60°,  source  depth  is  15  km. 

Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  24d:  Map  view  of  ray  bounce  points  froa  m  regional  model  that  includes  the 
Fielding  Moho  and  basin  databases.  Ray  take-off  dip  angle  is  60°,  source  depth  is  15  km. 
Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  25a:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  database.  Ray  take-off  dip  angle  is  70°,  source  depth  is  1  km. 

Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  25b:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  and  basin  databases.  Ray  take-off  dip  angle  is  70°,  source  depth  is  1  km. 
Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  26c:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  database.  Ray  take-off  dip  angle  is  70°,  source  depth  is  15  km. 

Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  25d:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  and  basin  databases.  Ray  take-off  dip  angle  is  70°,  source  depth  is  15  km 
Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  26a:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  database.  Ray  take-off  dip  angle  is  80°,  source  depth  is  1  km. 

Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  26b:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  and  basin  databases.  Ray  take-off  dip  angle  is  80°,  source  depth  is  1  km. 
Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  26c:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  database.  Ray  take  off-dip  angle  is  80°,  source  depth  is  15  km. 

Rays  were  shot  at  1°  azimuthal  increments  for  360°. 
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Figure  26d:  Map  view  of  ray  bounce  points  from  a  regional  model  that  includes  the 
Fielding  Moho  and  basin  databases.  Ray  take-off  dip  angle  is  80°,  source  depth  is  15  km 
Rays  were  shot  at  1°  azimuthal  increments  for  360°. 


3  Effects  of  Small  Scale  Structure 

3.1  Previous  work  and  limitations 

Ray  models  of  Lg  phase  are  useful  in  qualitatively  predicting  gross  features  of  blockage  by 
large  scale  crustal  transitions,  but  the  effects  of  smaller  scale  heterogeneity  are  more 
difficult  to  numerically  model.  The  most  important  effects  of  small  scale  heterogeneity  can 
still  be  accurately  included  in  modeling  by  the  Rayleigh-Born  approximation  to  scattering 
when  the  scale  length  of  the  heterogeneity  is  much  smaller  than  the  wavelength.  Realistic 
waveform  complexity,  however,  is  difficult  to  achieve  in  a  ray  method  without  fine  scale 
layering  and  large  numbers  of  scatterers.  The  coda  of  Lg,  including  equalization  of  Lg 
energy  on  all  components  of  motion,  has  been  attributed  to  back-scattering  of  seismic 
energy  off  lateral  heterogeneities  in  the  crust  (Aki  and  Chouet,  1975).  The  heterogeneities 
that  contribute  to  the  scattered  field  are  contained  within  two  ellipses  (figure  27a)  that  have 
foci  at  source  and  receiver  and  have  a  major  axes  equal  to  vt  and  v(t+dt),  where  v  is  the 
velocity  of  the  wave  field  and  t  is  time.  In  this  section,  the  effects  of  fine  scale 
heterogeneity  in  plane-layered,  laterally  homogeneous  structures,  are  modeled  by  the 
locked  mode  method  (Cormier  et  al.  1991,  Harvey  1981)  together  with  the  Rayleigh-Born 
approximation. 

The  Rayleigh-Bom  approximation  was  first  used  in  the  2-D  surface  wave  problem  by 
Kennett  (1972).  Snieder  (1986)  developed  the  theory  for  3-D  heterogeneities  by  using  a 
dyadic  representation  of  the  far-field  Green’s  function  in  a  laterally  homogeneous  medium. 
Wang  and  Hermann  (1988)  used  Born  scattering  to  investigate  surface  and  body  waves  in  a 
layered  medium.  Wang  and  Hermann’s  study  (1988)  accounted  for  mode  conversion, 
wave  type  conversion,  finite  scatterer  distribution  and  attenuation  from  scattering  as  well  as 
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from  intrinsic  anelasticity.  Wang  and  Hermann's  study  established  the  importance  of  mode 
conversion  and  found  that  Rayleigh  to  Love  wave  scattering  may  be  small  except  at  close 
distances.  Due  to  the  computational  enormity  of  the  problem,  their  study  limited  the  number 
of  modes  used  in  conversion  to  ten.  In  the  study  reported  here,  we  use  mode  summation 
similar  to  that  of  Wang  and  Hermann  (1988)  but  implement  Wu’s  derivation  for  Rayleigh 
scattering  in  terms  of  equivalent  sources  having  different  force  and  force  moments 
corresponding  to  perturbations  in  density  and  elastic  moduli.  This  enables  computation  of 
the  complete  mode  interaction  without  the  expense  of  computing  individual  mode  to  mode 
scattering  coefficients.  The  computation  of  modes  at  any  given  frequency  is  thereby 
reduced  from  m2  to  2m ,  where  m  is  the  number  of  modes,  thus  making  the  computation 
of  the  all  mode  interactions  a  tractable  problem  for  any  frequency  band. 

3.2  Modeling  method 

We  begin  by  following  Wu  and  Aki’s  (1985)  adaptation  of  the  Gubematis  et  al.  (1977a, 
1977b)  formulation  of  the  single-scattering  problem.  Given  an  arbitrary  heterogeneity  with 
variations  in  parameter’s 

A(x)  =  A0  +  8X(x),  <5A  «  A0 

p(x)  =  i u0  +  S/j.(x),  8p  « p0 

p(x)  =  p0  +  Sp(x),  where  8p  «  p0 

and  is  located  in  a  homogeneous  medium  with  parameters  X.0,  ji()  and  po.  The  displacement 
field  U  can  be  written  as  the  sum  of  the  direct  wave 

u=u°+u' 

Scatterers  are  treated  as  equivalent  sources  using  the  Bom  approximation  (Haddon  and 
Cleary,  1974;  Gubematis  et  al.,  1977  a,  b;  Aki  and  Richards,  1981;  Wu  and  Aki,  1985). 
For  Rayleigh  scattering,  in  which  the  wavelength  is  much  larger  than  the  inclusion,  the 
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scattered  field  can  be  equated  to  a  point  source  consisting  of  a  single  force  proportional  to 
the  density  contrast  and  a  moment  tensor  proportional  to  the  contrasts  in  elastic  moduli . 
The  equivalent  point  force  at  a  scattering  heterogeneity  can  be  written  as 


F^lQdV-tfSpVUffa) 

v 

and  the  moment  tensor  can  be  written  as 
=  j  Q,xtdV(x) 

V 

Mik = ~s,k  •  u°  (x0 )  -  m 


(3.1) 


r 


du?u») ,  auiM 


dxi 


dx: 


(3.2) 


where  Q.  is  the  equivalent  body  force,  V  is  the  volume  of  the  inclusion,  <5^  6fl,  and  $P  are 
the  average  perturbations  in  Lame’s  parameters  and  density,  and  8..  is  the  Kronecker  delta. 
It  may  be  noted  that  the  5)iV  in  equation  3.1  has  only  diagonal  components,  all  equal, 
which  corresponds  to  an  explosion.  The  moment  tensor  due  to  8(iV  is  a  symmetric  tensor 
that  may  have  both  diagonal  and  non-diagonal  elements.  The  diagonal  elements  of  5p.V 
correspond  to  on-line  force  couples,  and  the  off  diagonal  element  pairs  correspond  to 
torsion  double  couples. 


All  displacements  appearing  in  the  equivalent  source  representations  are  calculated  as  a  sum 
of  locked  Rayleigh  and  Love  modes,  using  a  modified  version  of  the  locked  mode  code  of 
Harvey  (1981).  The  displacement  at  the  inclusion  is  calculated  by  summing  over  modes  at 
each  frequency, 

n  mod  es 

U\k,z,(o)=  £ U°(k,z,co ) 

k=\ 

where  k  is  wavenumber,  z  is  depth  and  to  is  angular  frequency.  The  equivalent  point  force 
is  then  computed  for  each  inclusion  and  reradiated  towards  the  receiver.  This  entails  taking 
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the  horizontal  and  vertical  derivatives  of  the  stress  displacement  vector  to  compute  the  force 
moment  tensor  of  the  equivalent  point  force  vector.  The  scattered  far  field  (Aki  and 
Richards,  1981;  Wu,  1984)  is  then  calculated  as 

Ul  =  F  *  G  +  M  *G 

with  F  the  point  force,  M  the  equivalent  moment  and  G,  the  Green’s  function  response. 
This  convolution  is  performed  within  the  mode  summation  program,  by  looping  over  the 
modes  again  and  computing  tiie  response  due  to  the  point  force  and  moment  tensor.  The 
original  program  allows  for  sources  and  receivers  to  be  located  in  arbitrarily  in  a  Cartesian 
coordinate  system  and  rotates  this  coordinate  system  into  the  source-receiver  paths  so  that 
motion  can  be  partitioned  into  a  sum  of  either  pure  Rayleigh  and  or  pure  Love  modes. 
These  coordinate  rotations  are  performed  for  both  the  path  segment  from  the  source  to 
scatterer  and  the  path  segment  from  the  scatterer  to  receiver.  Figure  28  displays  a  map 
view  of  the  geometry  for  the  source-scatterer-receiver  used  in  computations. 

In  all  calculations,  a  step  function  time  dependence  is  assumed  for  the  displacement  in  the 
near-field,  corresponding  to  a  delta  function  time  dependence  for  both  P  and  S  body  waves 
observed  in  the  far-field.  This  is  assumed  to  be  a  good  approximation  for  earthquakes  of 
body  wave  magnitude  2.5  and  smaller,  having  far-field  P  corner  frequencies  higher  than  2 
Hz. 


3.3  Results  of  computational  experiments 

3.3.1  A  single  inclusion 

To  determine  where  in  the  wave  guide  it  is  necessary  to  know  the  detailed  fine  scale 
structure  that  affects  Lg  propagation,  a  simple  scenario  is  tested.  The  scattered  field  is 


66 


calculated  for  a  single  inclusion  for  a  source-receiver  pair  that  are  located  17.75  km  apart. 
The  structural  model  used  is  for  the  Pinion  Flats  region  in  California  (Scott  et  al.,  1994). 
This  region  is  chosen  due  to  a  well  constrained  velocity  model  from  prior  studies  as  well  as 
the  arrays’  proximity  to  a  large  regional  fault  that  may  have  scattering  effects.  The  structure 
from  the  region  used  to  calculate  the  mode  catalogue  is  in  table  2.  While  the  range  is  too 
close  for  a  true  Lg  phase  to  develop,  the  coda  following  the  S  arrival  that  is  the  result  from 
lateral  heterogeneities  has  the  same  characteristics  as  coda  Observed  at  greater  distances. 
The  scattering  bodies  are  assumed  buried  in  this  study.  This  simplifies  the  algorithm  and 
avoids  more  computations  of  perturbations  at  the  surface  boundary.  An  impulsive 
explosion  source  in  a  cavity  is  used  to  illustrate  the  conversion  of  modes  and  generation  of 
Lg.  An  explosive  source  in  a  homogeneous  medium  will  not  generate  energy  on  the 
transverse  component,  so  any  energy  observed  here  on  the  transverse  component  is  the 
result  of  mode  conversion  through  scattering  by  the  inclusion. 

The  volume  of  the  inclusion  used  in  this  study  is  5  km^  .  This  corresponds  to  a  sphere  of 
radius  of  approximately  1 .59  km.  For  a  frequency  of  3  Hz,  this  is  slightly  larger  than  the 
condition  that  provides  for  validity  of  the  Rayleigh  approximation: 


where  cmin  denotes  the  minimum  velocity  and  kmax  denotes  the  corresponding  wavenumber, 
0)  is  angular  frequency,  and  R  is  the  radius  of  the  volume.  Such  a  violation  in  the  volume 
size  may  produce  amplitude  errors  on  the  order  of  30%,  but  the  fundamental  scattering 
characteristics  will  valid.  The  variation  in  Lame’s  parameters,  |i  and  X,  used  for  the 
heterogeneity  is  -8%  and  density,  p,  is  -5%.  These  variations  are  used  to  represent  realistic 
variation  observed  in  crustal  seismic  velocities  and  densities. 


67 


The  inclusion  is  placed  at  a  constant  depth  of  2  km  and  then  moved  away  from  the  source 
in  2  km  increments  at  a  45  degree  azimuth.  The  geometry  for -this  is  illustrated  in  figure  29. 
The  inclusion  is  then  held  at  a  constant  range  of  2  km  and  45  degrees  and  it  is  placed  at 
depth  in  2  degree  increments.  The  same  procedure  is  used  to  test  heterogeneities  at 
distances  relative  to  the  receiver.  Synthetics  of  the  scattered  component  are  then  computed 
with  the  Bom-mode  method  and  plotted  for  each  location.  Figures  (30-33)  illustrate  the 
synthetics  that  resulted  from  this  test. 

The  synthetics  of  the  scattered  component  for  a  near-source  scatterer  at  constant  depth 
show  that  the  amplitude  in  the  group  velocity  window  2.75-3.5  km/s  rapidly  decreases 
with  distance  as  the  inclusion  is  moved  out  at  45  degrees  from  the  source  in  2  km 
increments  (figure  30).  The  transverse  energy  (figure  30b)  that  is  generated  is  on  the  same 
order  as  the  scattered  components  generated  on  the  radial  and  vertical  components  (figures 
30a,  c).  The  synthetics  generated  from  scatterers  at  a  constant  range  of  2  km  and  for  depths 
to  10  km  in  increments  of  2  km  are  plotted  in  figure  3 1 .  The  energy  decay  on  all 
components  appears  much  faster  with  depth  than  previously  observed  with  increasing 
range.  This  is  a  little  misleading  due  to  the  fact  that  the  source-scatterer-receiver  distance 
increases  more  rapidly  when  the  inclusion  is  moved  at  the  same  interval  in  the  z  direction  as 
in  the  horizontal  direction.  The  synthetic  seismogram  that  is  due  to  the  scattered  field 
illustrates  how  the  scattered  wave  field  amplitude  varies  with  inclusion  distance  and  depth 
from  the  source.  The  inclusions  closest  to  the  source  contribute  most  to  the  amplitude  of  the 
scattered  wave  field.  This  is  noted  in  both  the  horizontal  range  and  the  vertical  range. 

The  synthetics  from  models  where  the  heterogeneity  is  placed  near  the  receiver  yield  very 
similar  results.  In  both  the  increasing  horizontal  range  and  depth  cases  (figures  32-33), 
those  heterogeneities  closest  to  the  receiver  yield  the  highest  amplitude  scattered  wave.  It  is 
noted  that  transverse  energy  is  observed  for  all  placements  of  the  heterogeneity. 
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The  components  due  to  the  near-source  and  near-receiver  scattered  contributions  are 
summed  with  the  direct  arrivals  (figure  35a)  and  compared  to  the  direct  arrivals.  It  is  noted 
that  some  transverse  energy  is  generated,  but  not  to  the  same  extent  of  the  radial  or  vertical 
component.  The  vertical  and  radial  components  indicate  attenuation  due  to  deconstructive 
interference.  The  length  of  the  coda  or  wave  train  that  follows  is  lengthened  past  10 
seconds.  This  is  better  illustrated  in  the  envelope  of  the  traces  (figure  35b).  The  envelope 
function  is  defined  as 

env{t)  =  -Jx2(t)  +  y2(t) 

where  x(t)  is  the  original  signal  x(t)  and  y(t)  is  the  Hilbert  transform.  The  envelope  function 
is  a  useful  way  to  observe  decay  on  complex  signals  such  as  Lg  coda.  The  envelope  of  the 
direct  arrival  clearly  has  a  shorter  duration  and  decays  much  faster  than  the  envelopes  that 
include  all  of  the  scattered  components.  The  coda  in  this  case  was  significantly  increased 
with  the  addition  of  only  eighteen  point  scatterers,  thus  creating  a  more  complex  looking 
waveform.  In  this  case  the  majority  of  the  contribution  came  from  shallow  scatterers  (2  km) 
and  whose  range  was  increased  horizontally.  The  contribution  of  the  individual  inclusions 
rapidly  drops  with  increasing  depth,  as  can  be  seen  when  comparing  maximum  envelope 
amplitudes  for  the  Lg  group  velocity  window  3  3.5  km  /s  (figure  34). 


3.3.2  Multiple  inclusions 

The  crust  of  the  earth  is  thought  to  be  highly  heterogeneous  on  all  scales,  therefore,  in  this 
section,  clusters  of  inclusions  will  be  modeled  at  different  ranges  and  depths  to  determine 
to  what  extent  the  fine  scale  structure  needs  be  known  accurately  to  model  Lg.  Both 
shallow  earthquakes  and  explosions  (1  km  depth)  are  used  as  sources  to  generate  complete 
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synthetic  seismograms  and  to  examine  Lg  generation  in  an  effort  to  understand  what  effect 
the  source  type  has  on  scattering.  To  further  examine  range  and  depth  dependence, 
inclusions  are  put  arbitrarily  near  the  source,  near  the  receiver,  and  at  mid-range  spatial 
distributions.  Synthetics  are  generated  in  models  where  the  heterogeneities  are  at  shallow 
crustal  depths  of  3  and  5  km  and  deeper  crustal  depths  25  to  27  km  for  a  source-receiver 
distance  of  225  km  apart  (figure  36).  The  inclusions  have  variations  in  Lame’s  parameters 
(X  and  X  -9%  and  density,  o,  -3%,  respectively.  The  volume  of  the  inclusion  is  25  km3. 
The  parameters  are  chosen  on  the  upper  limit  of  what  is  expected  in  the  crust.  Scattering 
scales  linearly  with  the  volume  of  the  of  the  inclusions. 

The  first  source  type  considered  is  an  earthquake  source.  Source  parameters  are  a  step 
function  time  dependence  of  near-field  displacement  a  strike,  dip,  and  rake  of  61°,  70°, 

14°  respectively,  appropriate  for  an  earthquake  typical  for  the  Caucasus  region.  The  radial 
component  synthetics  due  to  the  earthquake  source  (figure  37)  show  the  contribution  from 
shallow  scatterers  (in  gray)  are  more  than  twice  the  size  of  the  contribution  from  deeper 
scatterers,  except  for  the  mid-ram  e  case.  Near-source  shallow  scattering,  however,  greatly 
exceeds  mid-range  shallow  or  deep  scattering  and  contributes  most  in  amplitude  to  the  final 
synthetic  (figure  38a-b).  The  duration  of  the  Lg  is  extended  as  scattered  energy  arrives  later 
in  the  Lg  group  velocity  window  when  compared  to  the  direct  arrival.  It  should  be  noted 
that  the  amplitudes  of  the  synthetics  due  to  scattering  are  plotted  relative  to  each  other  and 
not  relative  to  the  direct  wave,  as  the  parameters  of  the  inclusions  may  be  changed  to 
produce  greater  or  smaller  scattering  amplitudes  relative  to  the  direct  wave.  The  direct  wave 
is  plotted  for  shape  and  duration  reference  only.  In  the  case  of  the  transverse  component 
due  to  the  earthquake  source,  almost  all  of  the  scattered  energy  that  is  seen  is  due  to 
shallow  near-source  and  near-receiver  components  with  little  or  none  due  to  the  mid-range 
scattering  (figure  39).  The  scattered  energy  is  also  seen  to  be  peaking  later  in  the  Lg 
window  compared  to  the  Lg  in  the  direct  wave.  The  envelope  of  the  transverse  component 
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(figure  40b)  due  to  shallow  scattering  is  lengthened  considerably  and  extends  throughout 
the  Lg  group  velocity  window.  The  vertical  component  due  to-the  earthquake  source  does 
contain  considerable  energy  due  to  mid-range  scattering  (figure  41).  Near-source  shallow 
scattering,  however,  is  still  the  dominant  contribution  to  the  synthetic,  adding  most  of  the 
energy  later  in  the  Lg  group  velocity  window.  This  is  due  to  the  large  amplitude 
fundamental  mode  Rayleigh  wave  that  is  scattered  near  the  source  or  receiver  before  rapidly 
attenuating.  In  figure  42  a-b  the  direct  vertical  component  synthetic  is  plotted  along  with  the 
synthetics  of  the  combined  amplitudes  due  to  shallow  scattering  and  due  to  deep  scattering. 
The  envelope  of  the  vertical  synthetics  show  that  the  shallow  scattering  contributes  energy 
that  peaks  at  70  seconds  and  tails  off  near  the  end  of  the  group  velocity  window. 

The  same  scattering  geometry  is  used  to  compute  synthetic  seismograms  for  the  scattered 
energy  due  to  an  explosion  source  in  a  cavity  of .  1  km^.  The  radial  components  of  the 
direct  and  scattered  arrivals  are  plotted  in  figure  42-44.  The  shallow  near-source  scattering 
is  considerably  larger  than  the  other  synthetics.  Shallow,  near-receiver  scattering  also 
contributes  to  the  waveform,  but  inclusions  that  are  at  a  mid-range  distance  between  the 
source  and  receiver  produce  a  very  small  amplitude  scattered  wave  at  the  receiver  (figure 
44a).  The  same  test  for  deep  scatterers  produces  similar  results  for  range,  but  much  smaller 
than  from  shallow  scattering  (figure  44  b).  The  shallow  scattering  makes  up  almost  all  of 
the  scattered  energy  seen  when  summed  (figure  44  a,  b).  The  Lg  that  is  generated,  peaks 
around  74  seconds  and  decays  several  seconds  after  the  Lg  group  velocity  window.  The 
transverse  component  has  no  direct  wave  due  to  an  explosive  source.  The  transverse 
components  that  are  generated  from  scattering  are  plotted  in  figure  47  a.  It  is  obvious  that 
the  shallow  contributions  dominate  here,  as  well  as  the  near-source  and  near-receiver 
contributions.  The  synthetics  for  the  transverse  component  are  summed  and  plotted  (figure 
445  b).  The  Lg  that  is  generated  illustrates  realistic  Lg  shape  and  length.  The  coda  is  seen 
to  emerge  at  the  beginning  of  the  Lg  group  velocity  window  and  taper  off  a  few  seconds 
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after  it.  The  vertical  component  due  to  the  scattering  illustrates  the  same  general 
characteristics  (figure  46  a,  b).  There  is  a  significantly  larger  near-source  shallow  effects, 
with  shallow  near-receiver  scattering  contributing  as  well.  The  deep  near-source  scattering 
effect  does  appear  to  be  considerable.  It  is  small,  however,  relative  to  the  effects  seen  from 
the  shallow  scatterer.  The  summed  vertical  components  show  (figure  47  a,  b)  that  the 
shallow  scattering  again  dominates  the  shape  and  amplitude  of  the  scattered  synthetic.  The 
coda  shape  and  amplitude  are  extended  throughout  the  Lg  group  velocity  window  and 
decay  similar  to  observed  data.  It  is  noted  that  the  S  arrival  (approximately  at  50  seconds) 
has  gained  considerable  amplitude,  shape  and  complexity  of  its  coda. 


3.3.3  Scattering  at  CNET 

The  mode-Born  scattering  is  used  to  model  regional  data  from  the  Caucasus  array  (CNET) 
to  see  if  the  shallow  near-source  and  near-receiver  scattering  can  adequately  enhance  the 
direct  arrival.  The  wave  path  chosen  for  study  traverses  a  region  with  shallow  sediments 
and  that  is  highly  faulted.  This  path  does  not  undergo  any  significant  crustal  thinning  and 
thus  may  be  well  represented  by  a  layered  model.  The  chosen  path  is  from  an  event  on 
September  22,  1992  00:38:48.347,  located  at  42.94  N  and  44.06  E,  to  the  southwest  of  the 
array.  The  body  wave  magnitude  of  the  event  was  2.5 1  and  the  depth  was  reported  as  0 
km.  The  data  (figure  48)  shows  strong  Lg.  The  horizontal  components  are  rotated  to  the 
great  circle  path  and  put  into  radial  and  transverse  form  (figure  48  b).  The  data  is 
bandpassed  filtered  between  1  and  2  Hz  and  plotted  (figure  49  a).  An  envelope  of  the  traces 
is  taken  to  emphasize  the  arrivals  in  the  Lg  group  velocity  window  (49  b).  The  radial  and 
the  transverse  component  both  have  a  broad  Lg  arrival  towards  the  end  of  the  group 
velocity  window  that  has  a  high  frequency  content,  easily  seen  in  the  envelope  function. 


72 


The  synthetic  of  the  direct  arrival  generated  using  mode  summation  produces  a  waveform 
that  has  simple  peaks  when  the  envelope  of  the  trace  is  plotted  (figure  50  b).  The  arrivals 
for  the  P  and  the  Lg  are  approximately  the  same  as  seen  on  the  vertical  trace  of  the  data 
(figure  49  b),  indicating  that  the  layered  velocity  model  for  the  region  works  well  with  this 
particular  source-receiver  path. 

Scattering  is  now  added  to  try  to  match  the  complexity  seen  in  the  regional  wave:form.  The 
placement  of  the  near-source  scatterers  is  similar  to  that  used  in  (figure  37  a).  Only  shallow 
scattering  is  assumed  based  on  the  results  of  the  prior  studies.  Variations  in  Lame’s 
parameters  are  9%  with  a  variation  in  density  of  7%,  the  inclusion  is  assumed  to  be  a 
sphere  with  radius  0.9  km.  The  result  of  near-source  scattering  for  this  event  is  plotted  in 
figure  51.  The  synthetics  of  the  scattered  waves  from  shallow  near-source  scatterers  have 
significant  Lg  and  complexity.  The  envelope  function  is  plotted  to  emphasize  peaks  of 
energy  in  the  waveform.  The  peak  that  is  observed  in  the  later  part  of  the  Lg  group  velocity 
window  of  the  direct  synthetic  is  broader  and  is  similar  in  shape  to  the  observed  data.  The 
result  for  near-receiver  combined  with  near-source  scattering  does  not  noticeably  alter  the 
scattered  components  due  only  to  near-source  scattering. 

Arbitrarily  placing  scatterers  near  the  source  and  receiver  in  this  particular  case  dramatically 
improved  the  match  between  observed  data  and  synthetic.  While  there  is  little  doubt  that 
arbitrarily  placed  shallow  scattering  can  adequately  model  Lg  coda,  it  is  unknown  whether 
known  small  scale  crustal  features  may  be  incorporated  into  a  model  capable  of  producing 
realistic  coda. 
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3.3.4  Scattering  at  Pinion  Flat 


The  conventional  thought  is  that  seismic  coda  is  the  result  of  the  scattering  of  waves  from 
randomly  distributed  subsurface  heterogeneities  (Aki  and  Chouet,  1975;  Sato,  1977).  It  has 
been  suggested  that  regional  crustal  features  such  as  faults  may  be  a  source  of  scattering 
(Mitchell  1995),  thus  offering  a  deterministic  model  based  on  known  crustal  features 
instead  of  a  random  distribution.  Recently,  IRIS  broadband  data  from  Pinion  Flat 
Observatory  (PFO)  has  been  observed  to  have  off-azimuth  energy  that  was  not  attributed  to 
the  standard  back-scattering  from  random  lateral  heterogeneities,  but  attributed  to  the 
regional  north-west  trend  in  crustal  structure  (Wagner,  1996).  As  an  alternative  to  a  model 
consisting  of  randomly  placed  heterogeneities,  Wagner  (1997)  suggests  allowing  the 
heterogeneities  to  be  spatially  anisotropic.  Spatially  anisotropic  heterogeneities  cause 
spatially  anisotropic  scattering,  and  provide  a  characterization  of  subsurface  structure  that  is 
more  consistent  with  certain  geological  models  (Wagner  and  Langston,  1992; 

Wagner,  1996). 

To  test  a  scenario  of  spatially  anisotropic  heterogeneities,  mode-Bom  scattering  is  applied 
to  model  the  complexity  seen  in  the  regional  wave-forms.  A  fault-like  structure  is 
approximated  by  a  linear  zone  of  inclusions  surrounding  the  source  and  aligned  at  a  30 
degree  azimuth  with  respect  to  the  source  (figure  52).  The  inclusions  have  volumes  of  0.2 
km3  and  variations  in  Lame’s  parameters,  X  and  (X,  of  3%,  and  p  of  2%.  The  distance 
between  source  and  receiver  is  17.75  km.  The  direct  synthetic  due  to  an  explosive  source 
in  a  cavity  of .  1  km3  is  plotted  with  the  complete  waveform  that  includes  scattering  from 
the  ’fault’  (figure  53  a).  The  complexity  and  shape  of  the  synthetics  that  account  for  the 
spatially  anisotropic  scattering  provide  for  realistic  coda  in  shape  and  amplitude.  The 
envelope  function  of  the  direct  and  complete  waves  are  taken  to  illustrate  the  amplitude  of 
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the  coda  (figure  53  b).  It  is  noticed  that  the  explosion  generates  significant  transverse 
energy  when  scattered  from  the  fault.  The  coda  envelope  lengths  double  in  time  and  there  is 
slight  attenuation  in  peak  amplitude  of  the  vertical  component  due  to  deconstructive 
interference. 

To  determine  if  the  synthetics  could  be  analyzed  for  off  azimuth  energy,  a  simple  array  on 
the  same  order  of  Pinion  Flat  was  created  to  record  synthetics  (figure  54).  Frequency-wave 
number  (f-k)  analysis  was  performed  using  SAC  to  plot  azimuthal  characteristics  in 
arriving  energy.  The  time  windowed  S  coda  from  this  synthetic  array  is  shown  in  figure 
60  along  with  the  f-k  plot  (figure  55b).  The  maximum  radial  wavenumber  is  4  radian- 
seconds/km,  and  the  figure  shows  linear  contour  spacing.  There  is  an  identifiable  trend  of 
30°  azimuth  in  the  f-k  plot,  which  correlates  well  with  the  strike  of  the  fault . 

3.4  Summary  of  scattering  experiments 

Small  scale  heterogeneity  is  well  modeled  by  a  single  scattering  approximation  and  can 
account  for  the  duration,  complexity,  and  component  equalization  seen  in  Lg  coda.  The 
mode-Bom  method  of  calculating  single  scattering,  which  uses  an  equivalent  point  source 
to  reradiate  the  seismic  energy,  is  a  fast  and  efficient  way  to  predict  the  complexity 
observed  in  regional  seismic  coda 

For  shallow  sources  such  as  mining  blasts,  surface  faults  and  nuclear  explosions,  shallow 
scattering  statistics  close  to  the  source  and  receiver  are  shown  to  be  much  more  important  in 
controlling  Lg  complexity  than  detailed  knowledge  of  scattering  statistics  all  along  the 
crustal  waveguide.  This  supports  the  hypothesis  that  shallow  sediments  account  for  much 
of  the  lateral  scattering.  The  effect  of  shallow  small  scale  scattering  would  be  compounded 
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if,  as  in  the  first  part  of  this  report,  sediment  basins  trap  Lg  and  then  proceed  to  scatter  Lg 
by  small  scale  heterogeneities  within  the  basin. 


Scattering  near-source,  can  convert  the  strong  fundamental  Rayleigh  mode  to  higher 
modes.  The  fundamental  mode  Rayleigh  wave  is  known  to  attenuate  rapidly  and  travel 
relatively  small  distances,  and  thus  does  not  contribute  equally  throughout  the  scattering 
ellipse.  If  the  fundamental  mode  Rayleigh  wave  is  not  modeled  correctly,  or  if  shallow 
attenuation  (Q)  is  inaccurate,  the  scattering  from  the  fundamental  mode  will  also  be 
incorrectly  modeled. 


Explosive  sources  are  seen  to  generate  transverse  energy  due  to  near  source  and  near 
receiver  small  scale  heterogeneity.  This  is  a  particularly  important  point  in  terms  of 
discrimination,  suggesting  that  in  regions  of  thick  sediment,  with  perhaps  structural 
complexity  of  folds  and  thrusts,  a  small  explosion  may  generate  enough  transverse  energy 
to  be  appear  earthquake-like.  This  indicates  that  in  regions  with  complex  sediment 
structure,  particular  attention  must  be  paid  to  array  placement.  The  mode-Bom  method 
therefore  may  be  used  to  predict  and  optimize  array  locations  to  enhance  discrimination 
capabilities  in  a  given  area. 

The  size  of  the  heterogeneities  that  were  capable  of  adding  complexity  to  the  wave  form  and 
at  the  same  time  not  dominate  its  shape  were  on  the  order  of  .5  -  1 .0  km,  with  small 
variations  of  Lame’s  parameter’s  of  3%  well  within  maximum  crustal  variations  of  5-6%. 


Preliminary  tests  using  spatially  anisotropic  heterogeneities  work  well  in  regions 
characterized  by  a  dominant  trend  of  deformation.  With  a  general  knowledge  of  basin 
structure  as  well  as  trends  in  regional  geology,  it  may  be  possible  to  model  seismic  data 
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very  accurately  without  resorting  to  a  completely  statistical  descriptions  of  crustal 
heterogeneities. 

An  understanding  of  the  effects  of  small  scale  heterogeneities  and  their  variation  as  a 
function  of  large  scale  crustal  structure  and  tectonic  setting,  together  with  the  gross  effects 
of  phase  blockage  of  large  scale  structure,  will  form  an  important  knowledge  base  for  the 
transportability  of  seismic  methods  of  detection,  discrimination  and  yield  estimation. 
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Figure  27:  Scattering  ellipse,  with  foci  at  source  and  receiver  and  length  of  major  axis 
equal  to  vt  and  v(t+At),  where  v  is  velocity,  and  t  is  time.  Assuming  for  simplicity  that 
scattering  occurs  with  out  mode  conversion,  for  a  particular  mode,  the  shaded  region 
represents  the  region  that  contributes  to  a  particular  velocity  window.  In  the  example 
above,  the  Lg  group  velocity  is  highlighted.  The  shaded  elliptical  zone  is  the  region  that 
may  contain  heterogeneities  that  would  contribute  to  amplitudes  in  the  group  velocity 
window  of  3.0  to  3.5  km/sec. 
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Receiver 


/  01:  Source  -scatterer  relative  position 

/  02:  Scatterer-receiver  relative  position  used  in 

/  F2,F3. 

03:  Scatterer-receiver  relative  position  for  FI  at 
scatterer.  This  is  equal  to  it/2-02 

04:  Scatterer-receiver  relative  position  for  equiva¬ 
lent  moment  source.  This  is  equal  to  02  -  01 


Source 


Figure  28:  Geometry  of  source-scatterer-receiver  used  in  Mode-Born  program. 
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layer  no. 

depth  to  top  P-wave  velocit; 
(km)  |  (km/sec) 

S-wave  velocity 
(km/sec) 

density 

(g/cmA3) 

1 

0.00000 

5.53000 

3.23000 

2.40000 

2 

0.34500 

5.66000 

3.31000 

2.43000 

3 

1.34500 

5.79000 

3.39000 

2.46000 

4 

2.34500 

5.90000 

3.46000 

2.49000 

5 

3.34500 

5.97000 

3.50000 

2.52000 

6 

4.34500 

6.02000 

3.53000 

2.55000 

7 

5.34500 

6.05000 

3.54000 

2.58000 

8 

6.34500 

6.07000 

3.56000 

2.61000 

9 

7.34500 

6.10000 

3.57000 

2.64000 

10 

8.34500 

6.12000 

3.58000 

2.67000 

11 

9.34500 

6.14000 

3.60000 

2.70000 

12 

10.34500 

6.17000 

3.61000 

2.73000 

13 

11.34500 

6.19000 

3.62000 

2.76000 

14 

12.34500 

6.22000 

3.64000 

2.77000 

15 

13.34500 

6.25000 

3.66000 

2.78000 

16 

14.34500 

6.28000 

3.68000 

2.79000 

17 

15.34500 

6.31000 

3.69000 

2.80000 

18 

16.34500 

6.33000 

3.71000 

2.81000 

19 

17.34500 

6.36000 

3.73000 

2.82000 

20 

18.34500 

6.39000 

3.74000 

2.84000 

21 

19.34500 

6.42000 

3.76000 

2.85000 

22 

20.34500 

6.45000 

3.76000 

2.85000 

23 

21.34500 

6.47000 

3.79000 

2.87000 

24 

43.34500 

8.10000 

4.60000 

3.34000 

25 

69.34500 

8.60000 

4.90000 

3.50000 

26 

399.34500 

11.00000 

6.00000 

6.00000 

Table  2:  The  mode  catalogue  used  for  Pinon  Flats  synthetics  was  generated  using  this 
structure  model,  based  on  the  optimal  1-D  Anza  velocity  model  (Scott  et  al.,  1994). 


illustration  is  for  near  source  geometry;  the  near  receiver  geometry  is  set  up  exactly  the 


82 


5-01  St 


5-0 1  X 


6“0l  X 


—  in  cp  ir> 


6-01  X 


6-01  X 


6“0l  X 


6-01  X 


6-0 1  X 


6-01  X 


to  show  Lg  arrival  window. 


Depth  =  2  Km  -1  2  E  I  Depth  =  2  Km 


Depth  =10  Km  .  component;  B)  Transverse  component;  C)  Vertical 


Range 


Radial  component;  B)  Transverse  component;  C) 


86 


■  Near  source-increasing  depth 
Near  source-increasing  range 
■-EJ-  Near  receiver-increasing  depth 
--Q-  Near  receiver-increasing  range 


Figure  34:  Plotted  are  the  maximum  envelope  amplitudes  from  the  phase  velocity  group 
velocity  window  3-3.5  km/s  for  the  vertical  component  in  figures  5-8.  The  near  source 
and  near  receiver  scenarios  produce  amplitudes  that  vary  almost  identically  as  the  inclu¬ 
sion  is  moved  in  both  increasing  horiziontal  range  and  depth. 


X  IQ-7  X  ID-7  X  10-7 
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Figure  35:  a)  Reference  waveforms  plotted  with  summation  of  scattered  synthetics, 
b)  Envelope  of  reference  waveforms  and  scattered  components.  Notice  that  the  maximum 
amplitude  has  decresed  but  the  length  and  amplitude  of  coda  has  increased. 
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Figure  37:  Radial  component  of  mode  synthetics  due  to  an  earthquake  source.  This  illustrates  range  and  depth  effects  on  the  scattered 
components  (see  figure  39  for  location  geometry).  Top  trace  is  the  direct  arrival  and  is  plotted  to  scale  with  other  all  arrivals  except 
the  synthetic  due  to  shallow  near  source  inclusions.  The  shallow  near  source  effects  dominates  by  a  factor  of  10.  It  is  noted  that  mid 
source  scattering  contributes  on  the  same  order  to  the  Lg  phase  arrival  as  near  receiver  scattering. 
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Direct  Arrival 


Seconds 

I  Does  not  include  direct  arrival 


Figure  38:  Radial  component  of  mode  scattering  synthetics  due  to  an  earthquake  source, 
a)  Illustrates  1)  direct  waveform,  2)  summation  of  synthetics  from  shallow  scattering  all 
ranges,  3)  deep  scattering  all  ranges  and  4)  entire  geometrical  spectrum.  Shallow  scatter¬ 
ing  effects  dominate  in  the  synthetic  that  is  a  sum  of  all  contributions,  b)  A  plot  of  the 
envelopes  of  the  waveforms  in  figure  a.  Note  that  the  shallow  scattering  geometry  may 
increase  the  coda  in  amplitude  and  duration. 
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effects  on  the  scattered  components  (see  figure  39  for  location  geometry).  Top  trace  is  the  direct  arrival.  It  is  noted  that  shal¬ 
low  near  source  scattering  and  near  receiver  contributions  dominate  scattering  effects  and  that  scatterers  at  mid  range  geom¬ 
etry  do  not  have  a  strong  effect. 


92 


Figure  40:  Tranverse  component  from  Earthquake,  a)  1)  direct  waveform,  2)  summation 
of  synthetics  from  all  depths  and  ranges  ranges,  3)  deep  scattering  all  ranges  and  4) 
shallow  scattering  all  ranges.  Deep  scattering  effects  have  a  smoothing  effect  on  the  coda 
in  the  synthetic  that  is  a  sum  of  all  contributions,  b)  A  plot  of  the  envelopes  of  the  wave¬ 
forms  in  fig  a.  Note  that  the  scattering  adds  significant  amplitude  and  length  to  Lg  com¬ 
pared  to  the  direct  wave  and  that  while  deep  scattering  smooths  and  attenuates  the  Lg,  it 
does  not  change  the  characteristics  significantly. 
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near  source  and  shallow  near  receiver  scattering. 
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I  I  Does  not  include  direct  arrival  3  5  3  0 

Figure  42:  Vertical  component  mode  scattering  synthetics  from  an  earthquake  source,  a) 
The  1)  direct  arrival,  2)  direct  arrival  and  all  scattered  components,  3)  direct  arrival  and 
contribution  from  deep  scattering,  4)  direct  arrival  and  contribution  from  shallow  scatter¬ 
ing.  Note  that  the  shallow  scattering  is  significantly  larger  than  the  deep  scattering,  b) 

The  envelopes  of  the  synthetics  are  plotted  to  illustrate  the  lengthening  of  coda  due  to 
scattering  from  lateral  heterogeneities.  Shallow  effects  contribute  most  to  lengthening  the 
coda  envelope. 
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Figure  43  a,b:  Radial  component  of  mode  scattering  synthetics  due  to  an  explosion 
source,  a)  This  illustrates  range  effects  from  shallow  scattering  1)  direct  arrival,  2)  near¬ 
source,  3)mid-range,  4)  near-receiver,  b)  Range  effects  due  to  deep  heterogeneities, 
scattering  1)  direct  arrival,  2)  near-source,  3)mid-range,  4)  near-receiver.  The  shallow 
near  source  effects  dominates  by  a  factor  of  10.  It  is  noted  that  mid-range  scattering 
contributes  on  the  same  order  to  the  Lg  phase  arrival  as  near-receiver  scattering. 


Figure  44:  Radial  component  from  explosion  source,  a)  The  1)  direct  arrival,  2)  direct 
arrival  and  all  scattered  components,  3)  direct  arrival  and  contribution  from  deep  scatter¬ 
ing,  4)  direct  arrival  and  contribution  from  shallow  scattering.  Note  that  the  shallow 
scattering  is  significantly  larger  than  the  deep  scattering.  Shallow  scattering  effects 
dominate  in  the  synthetic  that  is  a  sum  of  all  contributions,  b)  A  plot  of  the  envelopes  of 
the  waveforms.  Note  that  the  envelope  due  to  the  shallow  scatterers  may  add  length  and 
amplitude  to  the  Lg  coda. 
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Seconds 

Figure  45:  Transverse  component  of  mode  synthetics  due  to  an  explosion  source,  a)  This 
illustrates  range  and  depth  effects  on  the  scattered  components  (see  figure  39  for  location 
geometry).  It  is  noted  that  shallow  near  source  scattering  and  near  receiver  contributions 
dominate  scattering  effects  and  that  mid-range  scatterers  are  not  seen  much,  b)  Envelope 
of  combined  scattering  contributions.  Envelope  shows  emergence  and  decay  typical  of 
that  observed  in  Lg. 
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,  :  “  Direct  Arrival  1-Direct  2-Near  source  3-Mid-range  4-Near  receiver 


Figure  46:  Vertical  component  of  mode  synthetics  due  to  an  explosion  source,  a)  This 
illustrates  range  effects  from  shallow  heterogeneities  on  the  scattered  components  (see 
figure  39  for  location  geometry),  b)  Range  effects  due  to  deep  heterogeneities.  The 
shallow  near  source  effects  dominates  by  a  factor  of  10. 
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direct  Arrival  f  Does  not  include  direct  arrival 

Figure  47:  Vertical  component  synthetics  from  earthquake  source,  a)  The  1)  direct 
arrival,  2)  direct  arrival  and  all  scattered  components,  3)  direct  arrival  and  contribution 
from  deep  scattering,  4)  direct  arrival  and  contribution  from  shallow  scattering.  Note  that 
the  shallow  scattering  is  of  same  order  as  the  deep  scattering,  and  most  of  the  deep 
scattering  is  from  near  source  effects,  b)  The  envelopes  of  the  synthetics  are  plotted  to 
illustrate  the  lengthening  of  coda  due  to  scattering  from  lateral  heterogeneities.  Shallow 
effects  contribute  most  in  lengthening  the  coda  envelope. 
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Figure  48:  a)  Raw  data  from  of  event  on  Sep  22,  1992  recorded  at  station  KIV  in  the 
Caucasus  network  b)  Horizontal  components  have  been  rotated  to  great  circle  path 
(radial  and  tranverse). 


.  v  Seconds 

b) 

Figure  49:  a)  Rotated  horizontal  components  and  band  passed  1-2  hz  data  from  of  event 
on  Sep  22,  1992  recorded  at  station  KIV  in  the  Caucasus  network,  b)  Envelope  function 
taken  of  the  rotated  and  filtered  data. 
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Figure  50:  a)  Synthetics  for  a  direct  arrival  for  the  Caucasus  event  on  September  22, 
1992.  Earthquake  source  at  a  depth  of  5  km  is  used  is  used,  b)  Envelope  function  of  the 
synthetics. 
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Figure  5 1 :  a)  Mode  syntheics  due  to  a  cluster  of  shallow  scatterers  near  the  source  for 
Caucasus  event  on  September  22,1992.  b)  Envelope  function.  Note  that  the  peak  in  the 
late  Lg  arrival  has  been  attenuated  and  more  closely  matches  real  data. 
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Figure  52:  a)  Mode  syntheics  due  to  cluster  of  shallow  scatterers  near  the  source  and 
receiver  for  Caucasus  event  on  September  22,1992.  b)  Envelope  function.  Note  that  the 
addition  of  more  scattering  (other  than  shallow  source)  does  not  noticebly  alter  the 
synthetics. 


105 


Figure  52:  Geometry  of  scatterers  used  to  calculate  synthetic  due  to  a  possible  fault  near  Pinon  Flats  Observa 
tory.  The  illustration  is  for  near  source  geometry,  the  near  receiver  geometry  is  set  up  exactly  the  same. 
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Figure  55:  a)  Synthetics  due  to  fault  geometry,  windowed  around  the  Lg  arrival  for  F-K 
analysis,  b)  F-K  plot  illustrating  that  energy  is  arriving  off  azimuth.  (Source-receiver 
azimuth  is  0°) 


4  Conclusions  and  Recommendations 

Ray  modeling  of  Lg  can  be  a  useful  reconnaissance  tool  for  predicting  the  efficiency  of  Lg 
propagation  in  regions  in  which  Moho  and  basin  topography  are  known  from  deep  seismic 
sounding  and  other  geophysical  observations  or,  at  the  minimum,  constrained  by  known 
surficial  geology.  Both  Moho  topography  and  basin  topography  can  be  important  in 
blocking  all  or  a  portion  of  the  Lg  coda,  with  basin  topography  primarily  affecting  shallow 
crustal  sources.  Predictions  of  Lg  efficiency  obtained  from  ray  modeling  are  useful  for 
design  of  array  siting  for  maximum  detection  capability  and  for  identifying  problematic 
paths  of  propagation  that  will  require  customized  regional  discriminants. 

The  complexity  of  regional  phases  is  strongly  controlled  by  scattering  by  small  scale 
structure  close  to  the  source  and  receiver,  including  the  equalization  of  Lg  energy  on  all 
components  of  motion,  regardless  of  source  type.  It  is  not  necessary  to  have  detailed 
knowledge  of  small  scale  heterogeneity  all  along  the  paths  of  regional  phases  to  be  able  to 
estimate  such  features  as  coda  shape  and  duration,  as  long  knowledge  of  heterogeneity 
statistics  are  known  close  to  the  source  and  receiver.  Anisotropy  of  heterogeneity  statistics 
can  be  important  to  include  in  regions  having  known  trends  of  deformation. 
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